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ABSTRACT 

Kepler has detected numerous exoplanet transits by precise measurements of stellar light in a single visible-wavelength band. 
In addition to detection, the precise photometry provides phase curves of exoplanets, which can be used to study the dynamic pro¬ 
cesses on these planets. However, the interpretation of these observations can be complicated by the fact that visible-wavelength 
phase curves can represent both thermal emission and scattering from the planets. Here we present a semi-analytical model 
framework that can be applied to study Kepler and future visible-wavelength phase curve observations of exoplanets. The model 
efficiently computes reflection and thermal emission components for both rocky and gaseous planets, considering both homo¬ 
geneous and inhomogeneous surfaces or atmospheres. We analyze the phase curves of the gaseous planet Kepler-7 b and the 
rocky planet Kepler-10 b using the model. In general, we find that a hot exoplanet’s visible-wavelength phase curve having a 
significant phase offset can usually be explained by two classes of solutions: one class requires a thermal hot spot shifted to one 
side of the substellar point, and the other class requires reflective clouds concentrated on the same side of the substellar point. The 
two solutions would require very different Bond albedos to fit the same phase curve; atmospheric circulation models or eclipse 
observations at longer wavelengths can effectively rule out one class of solutions, and thus pinpoint the albedo of the planet, 
allowing decomposition of the reflection and the thermal emission components in the phase curve. Particularly for Kepler-7 b, 
reflective clouds located on the west side of the substellar point can best explain its phase curve. We further derive that the 
reflectivity of the clear part of the atmosphere should be less than 7% and that of the cloudy part should be greater than 80%, and 
that the cloud boundary is located at 11 ± 3 degree to the west of the substellar point. Eor Kepler-10 b, the phase curve does not 
show a significant phase offset, and any model with a Bond albedo greater than 0.8 would provide an adequate fit. We suggest 
single-band photometry surveys could yield valuable information on exoplanet atmospheres and surfaces. 

Subject headings: radiative transfer — atmospheric effects — planetary systems — techniques: photometric 
— planets: Kepler-10 b — planets: Kepler-7 b 


1. INTRODUCTION 

A great number of exoplanets have been discovered by 
precise photometry. NASA’s Kepler spacecraft, monitor¬ 
ing 160,000 stars in the sky, discovered that more than half 
of the stars should host planets smaller than Neptune (e.g., 
Eressin et al. 2013; Howard 2013). The CHaracterizing ExO- 
Planet Satellite (CHEOPS), the Transit Exoplanet Sky Survey 
(TESS), and the PLAnetary Transits and Oscillations of stars 
(PLATO), designed to search for exoplanets around nearby 
bright stars using the same technique as Kepler, have recently 
been selected by ESA and NASA for launch within the next 
decade. The precise measurements of light curves of stars in a 
single visible-wavelength band will continue to be a dominant 
way to detect exoplanets, especially rocky exoplanets, in the 
coming years. 

Beyond planet detection, Kepler and future transiting ex¬ 
oplanet search missions may also provide valuable informa¬ 
tion on the nature of exoplanets via measuring the planets’ 
phase curves. The phase curves reveal longitudinal informa¬ 
tion regarding the planets’ atmosphere or surface. At visible 
wavelengths, a phase curve could illustrate how a planet re¬ 
flects stellar light and provide an effective way to study the 
condensed-phase particles in the planet’s atmosphere or the 
planet’s surface (e.g., Madhusudhan & Buitows 2012; Hu et 
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al. 2013; Heng & Demory 2013; Demory et al. 2013). The 
short-wavelength wing of thermal emission may extend to the 
visible wavelengths and affect the phase curve if the planet 
is hot enough (e.g., Rouan et al. 2011). It is therefore use¬ 
ful to learn whether a phase curve of combined reflection and 
thermal emission can place constraints on the atmosphere or 
surface of an exoplanet. 

The phase curve at visible wavelengths has been a powerful 
diagnostic tool to characterize Solar-System bodies, a method 
orthogonal to analyzing spectral features. The center-to-limb 
variation of the reflectivity continuum and major spectral fea¬ 
tures of Jupiter has yielded detailed information on the loca¬ 
tion and layering of its clouds (Sato & Hanson 1979). The 
phase curve of Venus alone has provided strong constraints 
on the refractive index of its cloud particles (Arking & Potter 
1968), which were later improved to effectively only allow 
spherical sulfate particles by measuring polarization (Hansen 
& Hovenier 1974). The photometric light curves of asteroids 
have been used extensively to study their shape and rotational 
properties (e.g., Torppa et al. 2003). 

Eor exoplanets, the phase curve characterization has been 
mostly limited to mid-infrared wavelengths due to greater 
planet-to-star flux ratios than at visible wavelengths (e.g., 
Knutson et al., 2007,2009a,b, 2012; Cowan et al. 2007,2012; 
Crossfleld et al. 2010; Lewis et al. 2013; Maxted et al. 2013; 
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oplanets show phase peaks before secondary eclipses, which 
have been interpreted as evidence for eastward displaced hot 
spots driven by super-rotating planet-encircling jet streams 
in the atmospheres (e.g., Showman et al. 2008, 2009). The 
visible-wavelength phase curves have been observed for one 
Jupiter-sized exoplanet before Kepler (Snellen et al. 2009), 
and recently for a number of exoplanets with Kepler data (e.g., 
Demory et al. 2013; Esteves et al. 2014). 

The interpretation of visible-wavelength phase curves with¬ 
out any spectral information, as is the case for Kepler observa¬ 
tions and also expected for future space photometry missions, 
are likely to be complicated by the following two factors: (1) 
a visible-wavelength phase curve that contains reflected light 
will depend on the spatial distribution of clouds on the planet; 
(2) a visible-wavelength phase curve may contain both re¬ 
flection and thermal emission from the planet, if the planet 
is highly irradiated. 

To interpret the visible-wavelength phase curves of exo¬ 
planets, we construct a simple semi-analytical model that can 
be applied to both Kepler and future visible-wavelength phase 
curve observations of exoplanets. Models at different levels 
of sophistication have been developed to study atmospheric 
circulation on exoplanets and interpret their thermal phase 
curves (e.g.. Showman & Guillot 2002; Showman et al. 2008, 
2009; Rauscher & Menou 2010,2012; Castan & Menou 2011; 
Wordsworth et al. 2011; Heng et al. 2011a, b; Cowan & 
Agol 2011; Heng & Kopparla 2012; Perna et al. 2012; Show¬ 
man & Kaspi 2013; Showman et al. 2013; Dobbs-Dixon 
& Agol 2013; Parmentier et al. 2013; Mayne et al. 2014; 
Rauscher & Kempton 2014; Kataria et al. 2014; Showman et 
al. 2014). However, because the Kepler observations contain 
only a single band, it is impractical to constrain the compo¬ 
sition of the planet’s atmosphere or surface. We approach 
the problem based on a physically motivated parameteriza¬ 
tion: instead of studying the detailed physical processes, we 
try to constrain several overarching physical parameters from 
the observations. This way we can compare planetary scenar¬ 
ios and shed light on general questions such as, whether the 
planet has an atmosphere, and whether the atmosphere has a 
standing circulation pattern and/or patchy clouds. 

This paper is organized as follows. We briefly outline pos¬ 
sible planetary scenarios important for phase curves in § 2, 
and then describe our semi-analytical model for interpreting 
the visible-wave length phase curve of exoplanets in § 3. In 
§ 4 we apply our interpretation tool to study the phase curve 
of the gaseous planet Kepler-7 b, and in § 5 we apply our in¬ 
terpretation tool to study the phase curve of the rocky planet 
Kepler-10 b. We discuss model degeneracies and suggest how 
they could be addressed with additional observations in § 6 
and conclude in § 7. 

2. PLANETARY SCENARIOS 

To make the model applicable for both gaseous and rocky 
exoplanets, we consider the following four planetary scenar¬ 
ios (see Figure [T] for a schematic illustration). These four 
scenarios are separated by whether a planet has an atmo¬ 
sphere, and whether a planet has a homogeneous reflecting 
layer. Whether a planet has an atmosphere is important be¬ 
cause the atmosphere could cause the greenhouse effect and 
raise the emitting temperature of the planet. Whether a planet 
has a homogeneous reflecting layer is also important because 
an inhomogeneous reflecting layer may cause a phase curve 
offset. 

The first and second scenarios are for bare-rock planets 


without atmospheres, and they differ by whether the surface 
is homogeneous. In the “homogenous surface” scenario, the 
planet is covered by a solid surface or a molten lava ocean. 
In the “inhomogeneous surface” scenario, the planet’s sur¬ 
face contains patches of differing albedo and/or thermal in¬ 
ertia. One example of the inhomogeneous surface scenarios 
is that only parts of the planet’s dayside are molten, forming 
a lava lake (e.g., Leger et al. 2011^ The molten surface and 
the solid surface may have different reflectivities. A key fea¬ 
ture that makes this scenario relevant is that the lava lake does 
not need to extend symmetrically with respect to the substel- 
lar point. This is an analog of a recent simulation on the ice 
coverage of tidally-locked ocean planets around M dwarf stars 
(Hu & Yang 2014). If a molten surface has a lower reflectiv¬ 
ity than a solid surface of the same material, an asymmetric 
lava lake may lead to a phase shift via reflection or thermal 
emission. 

The third and fourth scenarios are for planets with at¬ 
mospheres, and they differ by whether the atmosphere has 
patchy clouds. The distribution of the condensate particles, 
controlled by atmospheric circulation (e.g., Parmentier et al. 
2013), may be inhomogeneous, e.g., asymmetric with respect 
to the substellar point. The patchy cloud scenario is motivated 
by a similar scenario proposed to explain a post-occultation 
phase shift for Kepler-7 b (Demory et al. 2013). A cloudy 
patch of the atmosphere should have different albedo than a 
clear patch of the atmosphere - typically brighter - therefore, 
patchy clouds can induce phase modulation at visible wave¬ 
lengths. We do not separate planets having thick atmospheres 
and planets having thin atmospheres in the model, because we 
use a heat redistribution efficiency (see Section 3), rather than 
the surface pressure or the wind speed, as the primary model 
parameter. 

3. A GENERAL MODEL FOR ANALYZING VISIBLE-WAVELENGTH 
PHASE CURVES 

We aim to construct a general model to describe the phase 
curve of an exoplanet observed at visible wavelengths. The 
phase curve measures Tp/T/s as a function of the orbital phase, 
where Fp is the flux of the planet and Fs is the emission flux of 
the parent star. The planet’s flux comes from disk-integrated 
reflection and thermal emission, as 

Fp = Fp + Fp, (1) 

where Fp is the flux of stellar light reflected by the planet and 
Ft is the thermal emission flux of the planet. The reflection 
component depends on the albedo of the planet, and the ther¬ 
mal emission component depends on the efficiency of heat 
redistribution, and a potential greenhouse effect if the planet 
has an atmosphere. 

The model is designed for planets that have circular orbits. 
The planetary systems for which the model would be applied 
to would have their primary and secondary transits well mea¬ 
sured. The timing and duration of the transits constrain the or¬ 
bital eccentricity (Kallrath & Milone 1999; Barnes 2007; Ford 
et al. 2008; Kipping 2008; Demory et al. 2011). Therefore, 
one can confirm a planet to have a circular orbit before apply¬ 
ing the model. Most Kepler planets that have measurement of 
the secondary eclipse depths have circular orbits (Esteves et 
al. 2014). 

3.1. Reflection Component 

We model the reflection component in the following two 
ways. The first possibility is that the reflection component 
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Figure 1. Schematic illustration of the four planetary scenarios important for visible-wavelength phase curves. The figure shows the planet’s dayside, and the 
star mark indicates the location of the substellai' point. In the homogeneous surface scenario, the planet is covered by a solid surface or a molten lava ocean, 
and does not have an atmosphere. In the inhomogeneous surface scenario, the planet has a surface containing patches of differing albedo and/or thermal inertia, 
and does not have an atmosphere. One example of the inhomogeneous surface scenarios, as depicted on the figure, is that a fraction of the dayside is covered 
by molten lava (dark red) and the rest is covered by solid surface (light red). The solid part can be more reflective than the molten part. In the homogeneous 
atmosphere scenario, the planet has an atmosphere free of clouds or fully covered by clouds. In the patchy cloud scenaiio, a fraction of the dayside is cloud-free 
(dark grey) and the rest is cloudy (light grey). The cloudy part is more reflective than the cloud-free part. 


is symmetric with respect to the occultation. This cotTe- 
sponds to the “homogeneous surface” and “homogeneous at¬ 
mosphere” scenarios. The second possibility is that the reflec¬ 
tion component is asymmetric with respect to the occultation. 
This corresponds to a planet with different reflectivities at dif¬ 
ferent longitudes due to heterogeneous surface or atmosphere 
(i.e., the “inhomogeneous surface” and “patchy cloud" sce¬ 
narios). 

For symmetric reflection, we approximate the reflection 
component by the phase curve of a Lambertian sphere (i.e., 
all incident photons are isotropically scattered), scaled by a 
uniform geometric albedo (Ag) of the planet, namely 

FR(Ag) = Fs(—') Ag-[sin|a|-l-(7r-|a|)cos|a|], (2) 

Q, J TT 

where a is the phase angle of the planet that ranges from -tt 
to TT (a = 0 is at occultation and a = tt is at transit), Rp is the 
planet’s radius, and a is the semi-major axis. The Lambertian 
approximation is sufficient for current investigation because it 
is very difficult to extract a non-Lambertian component from 
a symmetric phase curve (Seager et al. 2000; Madhusudhan & 
Burrows 2012). If the planet has a homogeneously reflecting 
atmosphere, the shape of its reflection phase curve will closely 
resemble the Lambertian phase curve (Seager et al. 2000; Ca- 
hoy et al. 2010). For an airless rocky planet, we have used 
the Hapke planetary regolith reflection model to compute the 
shape of its phase curve based on the method described in Hu 
et al. (2012), and confirmed that the resulting phase curve will 
be sufficiently approximated by a Lambertian one, as long as 
the particle size of the regolith is less than 100 pm. 

For asymmetric reflection we consider some longitudes to 
have a high reflectivity and some longitudes to have a low re¬ 
flectivity. One could envision a lava lake picture - the molten 
surface has a lower reflectivity than the solid surface. One 
could also envision a “hole-in-a-cloud” picture - the atmo¬ 
sphere is clear and poorly reflective at some longitudes near 
the substellar point and cloudy and highly reflective at other 
longitudes. We denote the low reflectivity as ro and the high 
reflectivity as rQ + r\. Here a positive value for ri represents 
the increase in reflectivity caused by surface freezing or cloud 
formation. We further define [^ 1 ,^ 2 ] as the local longitud43 

* The local longitude is defined such that the substellar meridian is at the 
longitude of zero, the dawn terminator is at —7r/2, and the dusk terminator is 
at 7r/2. The “dawn” and “dusk” are defined for a planet of prograde rotation. 


range having low reflectivity, and then the local longitude 
ranges having high reflectivity are [-7r/2,^i] and [^ 277 r/ 2 ]. 
The reflection component is 

FR=Fi{2rQ/3) + Fs(^^ [ cos(a-0)cos((/>)r/^, 

\a J 3 TT 

(3) 

in which the first term is the Lambertian phase curve, char¬ 
acterized by ro, and the second term contains integration over 
the high-reflectivity longitudes visible to observers, character¬ 
ized by ri. The integration is expressed in the observer’s lon¬ 
gitude {(j>) such that the observer is at the direction of = 0. 
The relationship between the observer’s longitude and the lo¬ 
cal longitude is ^ = (p-a. The exact expression for the sec¬ 
ond term as a function of and ^2 is explicitly calculated and 
given in AppendixiAl 


3.2. Thermal Emission Component 

We calculate the thermal emission component in the Kepler 
band by multi-color blackbody emission based on a longitudi¬ 
nal distribution of temperature. The thermal emission photons 
can be either from the planet’s surface, or from a certain pres¬ 
sure level in the planet’s atmosphere, which we broadly refer 
to as the “photosphere”. The phase dependency of the planet’s 
thermal emission is computed by 


Ft = rI 


Bk[T{ a,9,9 co^p dOdp , (4) 


where Bk is the Planck function integrated over Kepler's 
bandpass, and the coordinates are specified by observer’s lati¬ 
tude and longitude (0, p), such that the observer is at the direc¬ 
tion of{9 = 0^(j) = 0). The thermal emission component at each 
phase angle a is controlled by the corresponding temperature 
distribution T{9,(j)). 

The temperature distribution is determined by interaction 
between irradiation, heat redistribution, and radiative cooling. 
Atmospheric circulation models have been developed for ir¬ 
radiated gas giants (e.g.. Showman & Guillot 2002; Show¬ 
man et al. 2009; Rauscher & Menou 2010, 2012; Heng et al. 
2011a, b; Perna et al. 2012; Showman et al. 2013; Dobbs- 
Dixon & Agol 2013; Mayne et al. 2014; Rauscher & Kemp- 
ton 2014; Showman et al. 2014) and super Earths (e.g., Castan 
& Menou 2011; Wordsworth et al. 2011; Heng & Kopparla 
2012; Kataria et al. 2014). Here we aim at fast calculation of 
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the temperature distribution that would enable parameter ex¬ 
ploration. Therefore, we adopt the semi-analytical model of 
Cowan & Agol (2011), in which the atmosphere is mimicked 
by a rigidly rotating “photosphere” subject to irradiation and 
radiative cooling The temperature distribution is written as 

T(a,9,cl))=fTo(0)V(e,O, ( 5 ) 

in which Tq is the temperature of the sub-stellar meridian, V 
is the thermal phase function that only depends on the local 
longitude and a heat redistribution efficiency (e), and / is a 
scaling factor to account for a possible greenhouse effect. V is 
computed by solving Equation (10) in Cowan & Agol (2011), 

^ = -(max(cosi,0)-'P‘^), (6) 

e 

and the sub-stellar temperature is 

7i) = 7i(^) ' (l-AB)‘/"cos(0)'/^ (7) 

where Ts and Rs are the stellar effective temperature and ra¬ 
dius, respectively, and is the Bond albedo of the planet. 

In this toy model, the heat redistribution efficiency is de¬ 
fined as the ratio between the radiative timescale and the ad- 
vective timescale, i.e.. 


^ — 'T'rad^adv? 


(8) 


where Trad is the radiative timescale of the photosphere, and 
the adveCtive frequency CUadv is CUadv = ^^photosphere te’orbit- 
cjphotosphere Worbit ^Tc the angular velocities of the photo¬ 
sphere and the bulk part of the planet in the inertial frame of 
reference, respectively. When |e| ^ 1 the longitudinal varia¬ 
tion of temperature will be small as heat redistribution is much 
more efficient than radiative cooling; when |e| <§; 1 the planet 
will be in local thermal equilibrium and the day-night contrast 
will be large. The sign of e, inherited from the sign of uj^dv, 
indicates the direction of the equatorial jets and the thermal 
phase shift: when e > 0 the photosphere is super-rotating with 
respect to the planet’s orbit, and the thermal phase shift is east¬ 
ward (i.e., the peak of the thermal phase curve appears prior 
to the occultation, as in the case of hot Jupiter HD 189733b; 
e.g., Knutson et al. 2007); when e < 0 the photosphere is 
sub-rotating with respect to the planet’s orbit, and the thermal 
phase shift is westward (i.e., the peak of the thermal phase 
curve appears after the occultation). If the photosphere can be 
treated as a single atmospheric layer, the radiative timescale 
would be 


^rad 


CpP 


( 9 ) 


where Cp is the heat capacity of the atmosphere, P is the 
pressure of the thermal emission photosphere, g is the grav¬ 
itational acceleration of the planet, and a is the Stephan- 
Boltzmann constant (Showman & Guillot 2002). 

Note that dynamical timescales other than the horizontal 
advection timescale may be important in controlling the tem¬ 
perature distribution on hot Jupiters; for example, in some 
cases, the vertical advection timescales and horizontal grav¬ 
ity wave propagation timescales exert a controlling influence 


^ As the thermal emission radiation in the Kepler band would almost en- 
tirely come from the equatorial region on the planet, a rotating photosphere 
is sufficient to describe the longitudinal variation of the photosphere temper¬ 
ature. 


on the day-night temperature pattern (Perez-Becker & Show¬ 
man 2013). In our above toy model, none of these dynamical 
processes are considered; Wadv should simply be viewed as a 
proxy for dynamical adjustment of the atmospheric tempera¬ 
ture structure, by whatever mechanism. 

The scaling factor / in Equation Q is included as a free 
parameter to mimic a possible greenhouse effect when there 
is an atmosphere. / = 1 corresponds to the homogeneous 
surface and inhomogeneous surface scenarios without green¬ 
house effects, and / > 1 corresponds to the planet that has an 
atmosphere with infrared-absorbing molecules. Such a scal¬ 
ing parameter is necessary because there could be infrared ab¬ 
sorbers in the atmosphere, such as CO and CO 2 , which do not 
contribute significantly to the opacity at visible wavelengths. 
In other words, Kepler could potentially probe a deeper, and 
presumably warmer, layer on the planet compared to the layer 
that reflects stellar light. Using a uniform scaling factor to ac¬ 
count for the possible greenhouse effect is of course coarse, 
as the radiative transfer processes in the atmosphere could be 
quite different from the dayside to the nightside (e.g.. Bur¬ 
rows et al. 2008). However, without knowing the details of 
the atmospheric composition, a scaling factor would be best 
suited for the purpose of this model. 

3.3. Linking Reflection and Emission Components 

At this point we can link the longitudinal variation of reflec¬ 
tivity to the longitudinal variation of temperature. We assume 
that the pattern for lava or clouds follows a longitudinal distri¬ 
bution controlled by the local temperature of the surface or the 
atmosphere (see Appendix iBl for justification). We introduce 
a single, physically motivated parameter called “conden¬ 
sation temperature” in the following, to describe the freezing 
temperature of the molten lava, or the condensation temper¬ 
ature of the condensable species in the atmosphere. When 
T({P{e,^) > Tc, the surface is molten, or the atmosphere is 
cloud free. When ToViefl) < the surface is solid, or the 
atmosphere is cloudy. ToV(efl) = Tc defines the longitudinal 
boundaries of the lava lake or the hole in the cloud, i.e., and 
^2- We drop the cos(0)^/"^ term in Tq (Equation|7]i to determine 
and ^ 2 , because any phase curve signal will be dominated 
by the atmospheric properties near the equatorial region that 
cos(6i) ^ 1. 

Putting these pieces together, we And that for the symmet¬ 
ric reflection scenarios (i.e., lava planet and homogeneous at¬ 
mosphere), the model is fully specified by three independent 
parameters: the Bond albedo of the planet (Ag), the heat re¬ 
distribution efficiency (e), and the greenhouse factor /. The 
geometric albedo and the Bond albedo are linked by the phase 
integral; and since we work with the Lambert sphere assump¬ 
tion in this paper, Ag = 2 AbI3. Eor the asymmetric reflec¬ 
tion scenarios (i.e., homogeneous surface and patchy cloud), 
the model needs two additional parameters: the condensation 
temperature (Tc), and a reflectivity boosting factor (k). The 
latter is defined as 

n = Kro. (10) 

Ab and k are two independent parameters, from which ro and 
ri can be calculated as 


fo ■■ 


Ab 

l + ^g'n' 


ri 


kAb 

1 


( 11 ) 
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where q' is the phase integral given in Appendix and Fig¬ 
ure [T] For physical plausibility, we verify each simulation to 
ensure ro-l-ri < 1, and reject all attempted models that do not 
satisfy this criterion. 

To summarize, our model computes the combined reflection 
and thermal emission from an irradiated exoplanet observed 
at visible wavelengths, based on three (assuming symmetric 
reflection) or five (assuming asymmetric reflection) parame¬ 
ters that describe physical processes on the planet. Our model 
covers a wide range of potential scenarios for exoplanets, with 
or without an atmosphere. 

3.4. Fitting to Light Curves 

The semi-analytical model can be used to fit observed light 
curves. Typically, we do not include the primary transit, as¬ 
suming that it has been used to derived key planetary and or¬ 
bital parameters, but we do include the secondary occultation 
in the fit. We use the model of Mandel & Agol (2002) for 
the shape of ingress and egress of the eclipse. This way, the 
phase curve characteristics including the occultation depth, 
the phase amplitude, and the phase offset can be directly de¬ 
rived from our fitting results. 

Since our model computes phase curves very fast, we can 
use the Markov-Chain Monte Carlo (MCMC) method to ex¬ 
plore the parameter space and determine the posterior pa¬ 
rameter distribution of the three or five parameters from the 
phase curve observations. In practice, we calculate 2 Markov 
chains, each containing 1 million steps, for each planetary 
scenario, with the MCMC method implemented as in Haario 
et al. (2006). This number of steps is sufficient for the con¬ 
vergence of the Markov chains for all parameters in the cases 
studied here, and we validate the convergence by compar¬ 
ing the two chains using the standard Gelman-Rubin statistics 
{R < 1.01 for all parameters; Gelman & Rubin 1992). The 
first half of each chain is considered the “burn-in” period and 
removed from the final results. The physically allowed ranges 
and the prior distributions of the model parameters are; Ab 
uniformly ranges in [0, 1]; e uniformly ranges in [- 00 , 00 ]; / 
uniformly ranges in [1,2]; 7], uniformly ranges in [200, 3000] 
K; and k uniformly ranges in [0,oo]. The actual ranges of 
parameters used in the fits to specific planetary scenarios may 
be narrower than these general ranges. 

Our model does not require knowledge of the stellar radius. 
With our formulation, the planetary phase curve {Fp/Fs) de¬ 
pends on Rp/Rs and a / Rs, both of which are uniquely derived 
from the primary transit (Seager & Mallen-Ornellas 2003). In 
practice, the precision for the measurements of Rp/Rs and 
a/Rs is on the order of 1% for Jupiter-sized planets whose 
transits and eclipses are detected by Kepler (e.g., Esteves et 
al. 2014), and the precision for the super-Earth-sized planet 
Kepler-10 b is better than ~ 3% (Batalha et al. 2011). It 
is therefore legitimate to not propagate the uncertainties in 
Rp/Rs and a/Rs to the fitted parameters, when the uncer¬ 
tainties of the fitted parameters are much greater than ^1%. 
The model also depends on the stellar effective temperature. 
The stellar temperature is not cancelled out because of non¬ 
linear dependency of the Planck function on the temperature. 
For well-characterized Kepler stars, such as Kepler-7 b and 
Kepler-10 b, the precision on the stellar temperature is well 
within 1%. We have tested our models and found 1% change 
in the stellar temperature would produce negligible change in 
the phase curve. 

In this paper we do not explicitly treat the doppler beaming 
or ellipsoidal effects to the phase curve, assuming that these 
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Figure 2. The best-fit models for the phase curve of Kepler-7 b. The ma¬ 
genta lines are the modeled phase curves, and the other lines show contribu¬ 
tion of thermal emission (red), symmetric reflection (green), and asymmetric 
reflection (blue). The upper-panel shows the best-fit model for a planet hav¬ 
ing a homogeneous atmosphere, and the modeled phase curve has significant 
contribution from a thermal emission component characterized by a hot spot 
shifted westward. By definition in this scenario the asymmetric reflection 
component is zero. The lower-panel shows the best-fit model for a planet 
having an atmosphere with patchy clouds, and the modeled phase curve has 
significant contribution from the asymmetric reflection component attributed 
to the clouds. The clouds concentrate on the west side of the substellar point, 
driven by a hot spot shifted eastward. 

effects can be effectively removed for transiting planets be¬ 
cause relevant orbital parameters are known from the transits. 
The interaction between these effects and the atmospheric sig¬ 
natures will be discussed in a separate paper (Shporer & Hu 
2015). 

4. APPLICATIONS TO KEPLER-7 B 

We now apply our model framework to analyze the visible- 
wavelength phase curves of exoplanets. Kepler has provided 
data to derive visible-wavelength phase curves for a number 
of Jupiter-sized exoplanets, and among these planets, Kepler- 
7 b has the best signal-to-noise ratios for the occultation depth 
and the phase offset (Demory et al. 2013; Esteves et al. 2014). 
Kepler-7 b has a circular orbit with an orbital eccentricity less 
than 0.02, the 3 - cr upper limit derived from the light curve 
(Demory et al. 2011). As an example, we apply our model to 
analyze the phase curve of Kepler-7 b. We use the phase curve 
data of Demory et al. (2013) in this study. To focus on testing 
what we could learn from visible-wavelength phase curves, 
we only use the phase curve for the fit, and then consider the 
constraints of the occultation depths measured at longer wave¬ 
lengths by Demory et al. (2013). We show our results in Table 
[Hand Figures 1212 and|4] 

A bimodal distribution for the htted and derived parameters 
emerge when we apply the model framework to explain the 
phase curve of Kepler-7 b. The key feature of the observed 
phase curve is that the peak of planetary light occurs after the 








6 


Hu et al. 


Table 1 

Estimation of parameters for the gaseous planet Kepler-7 b, based on fitting to the observed phase curve. We use the photometry derived by Demory et al. 
(2013). In the general fit, all five parameters are allo\ved to vary in their physically plausible ranges. In the homogeneous atmosphere fit, the asymmetric 
reflection component is not included in the calculation, and then the model no longer depends on the cloud condensation temperature or the cloud reflection 
boosting factor. In the patchy cloud fit, we set e > 0, assuming the planet to have super-rotating equatorial winds. Under such assumption, patchy clouds are 

required to explain the observed phase curve (see the text). 


Parameter 

General 

Homogeneous Atmosphere 

Patchy Cloud 

Fitted Parameters 

e 

f 

T(K) 

K, 

-3.3 ~ 62 

1 1 q - K ).05 
^•^^- 0.16 

1476 2604 

14-. 82 

0.18 + 0.03 

-3.2 + 0.5 

1.23 + 0.02 

0.42 + 0.01 

20~82 

1480+10 

13~69 

Derived Parameters 
Occultation Depth (ppm) 
Phase Amplitude (ppm) 
Phase Offset (degree) 

^0 

^1 

^1 (degree) 

^2 (degree) 

41.5);” 

48.0±2.1 

'in 9 + 2.6 
- 3.0 

0.003 ~ 0.049 
0.15 ~ 0.94 
-13.8 ~0 

0 90 

44.5 + 2.0 

49.4+1.4 

3&.Q+1* 

39.0flj 

46.6+1.3 

36.2 + 2.7 

0.014 ~ 0.072 
0.92 + 0.04 
-11.2 + 2.7 

90 

Fit Quality 

Minimum /do f 

0.992 

0.991 

1.028 

BIC 

1340.8 

1327.1 

1387.1 


secondary occultation, i.e., a post-eclipse phase offset. Our 
model can generate this offset and provide a satisfactory fit to 
the observation (Figure |2]i. When all parameters are allowed 
to vary in their physically plausible ranges (i.e., the general 
fit), the posterior distribution shows two clearly separate pa¬ 
rameter spaces that produce the fit to the observation (Figure 
[2and|4|i. The bimodal distribution is especially apparent for 
the Bond albedo, which can take values around either 0.2 or 
0.4. Such a bimodal distribution indicates that two classes of 
models can be consistent with the observed phase curve. 

The two classes of models correspond to the homogeneous 
atmosphere scenario and the patchy cloud scenario, respec¬ 
tively. To separate the two classes, we additionally perform 
two fits, each corresponding to one planetary scenario, by lim¬ 
iting the ranges in which the fitted parameters can vary. For 
the homogeneous atmosphere fit, we assume the asymmetric 
reflection component to be zero, which makes the cloud con¬ 
densation temperature and the cloud reflection boosting factor 
dummy parameters. For the patchy cloud fit, we assume the 
planet to have super-rotating equatorial winds (e > 0). Un¬ 
der such assumption, homogenous atmosphere can no longer 
provide a fit, and patchy cloud would have to be invoked to ex¬ 
plain the post-eclipse phase offset. Table [T] shows the results 
of the two additional fits, and Figure [3 compares the poste¬ 
rior distributions from the two separate fits with those from 
the general fit, and shows that the two separate fits correctly 
capture the two classes of solutions suggested by the general 
fit. 

While producing the same phase offset, the two classes of 
models result in slightly different occultation depths, because 
they differ in the shape of the phase curve. Table [T] shows 
that the derived occultation depths differ by 2-a between the 
homogenous atmosphere fit and the patchy cloud fit, and Fig¬ 
ure [3 shows that the general fit would give the combination 
of the above two fits, and therefore have much greater error 
bars. Comparing with Demory et al. (2013), the occultation 
depth from the patchy cloud fit is 1-cr consistent, and that from 
the homogeneous atmosphere fit is only 2-a consistent. Sim¬ 


ilar bimodal distribution is also found for the posterior of the 
phase amplitude, but the two peaks differ to a lesser extent. 
This difference identified here highlights the importance to 
use appropriate models for the phase curve, even when the 
occultation depth bears most interest. A detailed look at the 
best-fit model phase curve (Figure |2| would reveal that the 
homogeneous atmosphere model has more thermal emission 
contribution than the patchy cloud model, and then the onset 
of the planetary light is smoother as the planet rotates from the 
nightside to the dayside. To fit the phase curve that contains 
the secondary occultation, the model automatically adjust the 
“zero” point of the planetary flux to seek the minimum 
which affects the determination of the occultation depth. 

4.1. Homogeneous Atmosphere 

The best-fit model for the homogeneous atmosphere sce¬ 
nario has thermal emission as the dominant source of plan¬ 
etary light. In this scenario, the post-occultation phase shift 
can be explained by a hot spot located on the west side of the 
substellar point (Figure |2]l. All three fitting parameters (the 
Bond albedo, the heat redistribution efficiency, and the green¬ 
house factor) are tightly constrained in this scenario (see the 
red lines in Figure O. With only three parameters, the sce¬ 
nario achieves a superior goodness of fit indicated by and 
significantly better value for the Bayesian Information Crite¬ 
rion (BIC) than any other scenarios (TabIe[T]). 

Assuming a homogeneous atmosphere, the visible- 
wavelength phase curve contains enough information to de¬ 
termine the three fitting parameters. The phase curve depends 
on these parameters non-linearly, and these parameters appear 
to be correlated but not fully degenerate (Figure|4|. The most 
prominent correlation is between the greenhouse factor (/) 
and the Bond albedo (A^), because the modeled temperature 
distribution is proportional to /(I -A^)'/'*. These two parame¬ 
ters tend to be correlated, in order to maintain the temperature 
and then the thermal emission component (Figure]?]). We also 
And that for a greater Ab, the model requires a greater (i.e.. 
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Condensation Temperature 7^ Cloud Reflection Boosting Factor k 

Figure 3. Posterior probability disti'ibution of the parameters for Kepler-7 b from MCMC simulations to fit the phase curve. Results from the general fit, the 
homogeneous atmosphere fit, and the patchy cloud fit are shown by different colors, and the color designation is tabulated on the lower right. The posterior 
probability distribution resulted from the general fit shows bimodal solutions, and appears to be the combination of the results from the homogeneous atmosphere 
fit and the patchy cloud fit. The homogeneous atmosphere fit tightly constrains the heat redistribution efficiency to —3.2 ±0.5, and the greenhouse factor to 
1.23 ± 0.02. The patchy cloud fit allows a positive e, but the phase curve yields no constraint on the exact value for e; instead, the phase curve well constrains the 
cloud condensation temperature. 


more negative) value for the redistribution efficiency e (Fig¬ 
ure Si. This is because when Ag is greater, the symmetric re¬ 
flection component becomes more significant and the thermal 
emission component becomes less significant. To keep the 
phase offset consistent with the observation, the offset of the 
thermal emission component alone would have to be greater, 
and therefore |e| would have to be greater. Furthermore, when 
|e| is greater, the temperature becomes lower at the day side 
and higher at the nightside, which in turn would affect the oc- 
cultation depth and the phase amplitude. In all, three “model- 
independent” observed quantities, the occultation depth, the 
phase amplitude, and the phase offset can uniquely determine 
the three fitting parameters. Once their values are found, 
we can separate the contribution from thermal emission and 
that from reflection, solely from the visible-wavelength phase 
curve. 

Particularly for Kepler-7 b, we find that the heat redistri¬ 
bution efficiency (e) is smaller than zero by 6 - cr, implying 
that the advective frequency, Wadv = <^photosphere “ ^orbit must 
have a large negative value, in order to explain the post- 
occultation phase shift. Therefore, the westward offset of the 
bright spot would seem to suggest that the air westward of 


the substellar point has the greatest temperature, and there¬ 
fore that air at photospheric pressures travels westward in a 
synchronously rotating reference frame. This could occur 
on a synchronously rotating planet if the photospheric-level 
winds were westward at low latitudes. However, to date, cir¬ 
culation models of highly irradiated, tidally locked exoplan¬ 
ets have generally predicted eastward equatorial winds and 
therefore eastward offsets of hot spots relative to the substel¬ 
lar point (e.g.. Showman & Guillot 2002, Cooper & Show¬ 
man 2005, Showman et al. 2008, 2009; Rauscher & Menou 
2010, 2012; Heng et al. 2011; Perna et al. 2012). Such east¬ 
ward jets were explained in a theory presented by Showman 
& Polvani (2011), which shows that the day-night thermal 
forcing induces global-scale waves that transport prograde an¬ 
gular momentum to the equator, allowing such a so-called 
“superrotating” equatorial jet to emerge. To date, no models 
of highly irradiated, synchronously locked planets have been 
published that exhibit a strong westward jet at the equator, as 
would seem to be needed to explain the westward offset in 
the Kepler-7 b phase curve in the homogeneous atmosphere 
scenario. 

In principle, non-synchronous rotation could contribute to 
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Figure 4. Correlations between fitted parameters for Kepler-7 b. The results from the general fit are shown as black dots. The figures show two separate 
populations con'esponding to the homogeneous atmosphere solution and the patchy cloud solution. Small Bond albedos {Ag), negative heat redistribution 
efficiencies (e), and greater-than-unity greenhouse factors (/) correspond to the homogeneous atmosphere solution; the other two parameters are unconstrained. 
The three constrained parameters are correlated, and their congelations are highlighted by inserted panel 1 and 2 that provide zoom-in views. Large Bond albedos, 
positive heat redistribution efficiencies, well-constrained condensation temperatures {Tc) near 1500 K, and positive reflectivity boosting factors (ac) correspond to 
the patchy cloud solution. The Bond albedo and the condensation temperature are correlated and their correlations are highlighted by inserted panel 3. 


a hot spot offset. Circulation models of non-synchronously 
rotating hot Jupiters have been computed by Showman et al. 
(2009), Showman et al. (2014), and Rauscher & Kempton 
(2014), considering rotation that is prograde (i.e., in the same 
direction as the orbital motion), with rotation periods both 
shorter and longer than the orbital period. Like their syn¬ 
chronously rotating counterparts, almost all of these mod¬ 
els develop fast eastward equatorial jets - fast enough to 
cause eastward motion in the synchronously rotating refer¬ 
ence frame regardless of the rotation rate - and therefore east¬ 
ward hotspot offsets. That said, one slowly rotating non- 
synchronous simulation in Rauscher & Kempton (2014) de¬ 
velops a “westward” equatorial jet that causes the thermal 
hotspot to be shifted west of the substellar point. In princi¬ 
ple, planetary rotation that is retrograde (i.e., in the direction 
opposite the orbital motion) could also lead to thermal hotspot 
offsets of the correct sign to explain the Kepler-7 b data. How¬ 
ever, given the short tidal spindown timescales for very hot 
Jupiters like Kepler-7 b, it is likely that the planet is close to 
synchronous rotation, in which case one would expect the ex¬ 
istence of a superrotating jet and an eastward hotspot offset. 

For Kepler-7 b, the dayside equilibrium temperature for the 


best-fit homogeneous atmosphere model is 1820 K, which is 
higher than the 3-cr upper limit of the brightness temperature 
measured at the Spitzer 3.6 pm band (Demory et al. 2013). 
With the caveat that Kepler and Spitzer may probe different 
pressure levels and have different brightness temperatures, the 
homogeneous atmosphere scenario appears to be inconsistent 
with the Spitzer observations. Based on the atmospheric cir¬ 
culation model results and the Spitzer observations, the ho¬ 
mogeneous atmosphere scenario is unlikely, which effectively 
makes the patchy cloud scenario to only plausible scenario for 
Kepler-7 b. 

4.2. Patchy Cloud 

The best-fit model for the patchy cloud scenario has the 
asymmetric reflection component as the dominant source of 
planetary light (Figure |2]). In this scenario, reflective clouds 
located on the west side of the sub-stellar median can best ex¬ 
plain the post-eclipse phase shift of Kepler-7 b. This is con¬ 
sistent with the explanation proposed by Demory et al. (2013) 
but our analysis offers more information on the atmosphere’s 
properties. The Bond albedo is well constrained (Table [T]), 
and is much greater than that in the homogeneous atmosphere 
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scenario. With a Bond albedo of ~ 0.4, the dayside equilib¬ 
rium temperature is consistent with the Spitzer observations. 

The asymmetric reflection component is produced by a 
cloud distribution, in which the east side of the sub-stellar 
median is devoid of reflective clouds. To form such a cloud 
distribution, the model requires non-zero positive values for 
the heat redistribution efficiency (e) and the cloud reflectivity 
boosting factor (k) (Table [Til. With a positive value for e, the 
east side is hotter than the west side; and if the cloud conden¬ 
sation temperature is suitable, condensation can only occur 
on the west side but not on the east side. This forms a “hole” 
in the cloud on the east side. Particularly for Kepler-7 b, the 
cloud boundary would have to be located at ~ 10 ° to the west 
of the substellar point to produce the observed phase offset 
(Table [T]i. A positive value for e indicates super-rotating jets 
that transport heat towards east, consistent with atmospheric 
circulation theories (Showman & Guillot 2002; Showman et 
al. 2011 ). 

Also, to explain the phase offset, the cloudy part of the at¬ 
mosphere must be more reflective than the cloud-free part, 
and our analysis shows that this reflectivity contrast would 
have to be quite significant to explain the phase curve of 
Kepler-7 b. We And that the cloud-free part of the atmosphere 
must be quite dark, having a reflectivity less than 7%, and the 
cloudy part of the atmosphere must be quite bright, having a 
reflectivity greater than 90%, at visible wavelengths (Table[T]i. 

The significant contrast in reflectivities is one of the main 
reasons why the model can constrain the Bond albedo. The 
Bond albedo cannot be too small, otherwise there would not 
be enough planetary light to explain the phase amplitude. Al¬ 
ternatively, the phase amplitude could be explained by ad¬ 
ditional thermal emission (i.e., by increasing the greenhouse 
factor /), but that would drive the overall phase shift towards 
the opposite direction than the asymmetric reflection. The 
Bond albedo cannot be too large also, because the hole in 
cloud has to be large enough (i.e., covering a significantly 
large part of the dayside), and dark enough for the signifi¬ 
cant contrast between the cloudy part and the cloud-free part. 
In all, the phase amplitude and the phase offset together put a 
tight constraint on the Bond albedo of the planet in the patchy 
cloud scenario. 

In addition to the Bond albedo, fitting to the Kepler phase 
curve constrains the cloud condensation temperature and may 
imply the physical properties of the condensate species. The 
cloud condensation temperature is correlated with the Bond 
albedo, as Tc oc (l-A^)'/^ (Figure |4]i. This is because the 
cloud condensation temperature needs to have a value be¬ 
tween the maximum and the minimum temperature on the 
dayside, in order to produce a patchy cloud distribution. Cau¬ 
tion should be exercised when comparing Tc with the conden¬ 
sation curves of potential cloud-forming materials. We com¬ 
pare Tc with the equilibrium temperatures in the model setup. 
However, the cloud may be located deep in the atmosphere 
at the pressure of 0 . 1-1 bars, and the true cloud condensation 
temperature may be higher than Tc- A more realistic range for 
the cloud condensation temperature would be between Tc and 
Tcf, where / is the derived greenhouse factor. For Kepler- 
7 b, we And this range corresponds to 1480-1730 K, using 
the values tabulated in Table [T] The condensation curves of 
Fe, Mg 2 Si 04 , MgSi 03 , and Cr cross the inferred temperature 
range at the pressure of 10 “^ -1 bars for the solar abundances 
(Fodders & Fegley 2006). Fe and Cr are strongly absorp¬ 
tive in the visible wavelengths, and thus cannot lead to the 


required high reflectivities. Mg 2 Si 04 and MgSi 03 are highly 
reflective (Sudarsky et al. 2003), and therefore are candidate 
cloud-forming materials for the atmosphere of Kepler-7 b. 

Finally, the efficiency of heat redistribution cannot be suf¬ 
ficiently constrained by the phase curve. For Kepler-7 b, a 
non-zero positive value for e is required, but thephase curve 
does not prefer a specific value for e (Figure |3|l - meaning 
that the phase curve is rather insensitive to the exact temper¬ 
ature distribution as long as a hot spot offset exists. This is 
because when e is large, the longitudinal variation of temper¬ 
ature would be small (Cowan & Agol 2011), and a slight ad¬ 
justment of Tc would be enough to keep the cloud boundary 
unchanged. This is different from the homogeneous atmo¬ 
sphere scenario, in which the phase offset directly depends on 
the hot spot offset and e can be tightly constrained. In the 
patchy cloud scenario, due to the uncertainties of the cloud 
condensation temperature, the phase curve sets a 1 - cr lower 
bound of 20 and does not yield an upper bound. In other 
words, the phase curve only requires a “fairly significant” heat 
redistribution, but cannot yield quatitative constraints on this 
parameter. 

5. APPLICATIONS TO KEPLER-10 B 

We derive a full phase curve of Kepler- 10b based on Kepler 
observations during the quarters 1 to 17. Kepler-10 b is a 1.4- 
R 0 , 4 . 6 -M 0 predominantly rocky planet (Batalha et al. 2011; 
Fogtmann-Schulz et al. 2014), making our result among the 
first revelation of any rocky exoplanets’ phase curve signa¬ 
tures. This observation is made possible by continuous mon¬ 
itoring of the system by Kepler that brings down the error 
budget of photometry. We also benefit from the fact that the 
planet is hot enough to have significant thermal emission con¬ 
tribution to the Kepler band (Rouan et al. 2011), and the star 
is intrinsically quiet (Batalha et al. 2011). We apply our semi- 
analytical model to analyze the phase curve of Kepler-10 b. 

5.1. Data Reduction 

Our data reduction is similar to the one for Kepler-7 b pre¬ 
sented in Demory et al. (2013). We use Kepler (Batalha et 
al. 2013) long-cadence simple aperture photometry (Jenkins 
et al. 2010) obtained during the quarters 1 to 17. We take 
into account the crowding matrix correction factor indicated 
in each FITS file on a quarter-per-quarter basis. We mitigate 
instrumental systematics by fitting the first four cotrending ba¬ 
sis vectors (CBV) to each quarter using the PyKE software 
(Still & Barclay 2012). We then normalize each quarter to the 
median. We account for photometric trends longer than four 
times the planetary orbital period by fitting a second-order 
polynomial to the out-of-eclipse data in the MCMC frame¬ 
work detailed below. We estimate and include the corrected 
noise the same way as in Demory et al. (2013). We And a 
nominal level (less than 10 %) of correlated noise throughout 
the dataset. 

5.2. Model-Independent Analysis 

Before isolating the planetary phase curve signal we search 
for all frequencies in the dataset to assess any risk of contam¬ 
ination. A typical Fomb-Scargle periodogram is not optimal 
in the case of datasets spanning long observations as slightly 
changing periodicities damp amplitudes in the power spec¬ 
trum. To quantify how frequencies and amplitudes evolve in 
our dataset, we perform a wavelet transform analysis (e.g., 
Torrence & Compo 1998) using the weighted wavelet Z- 
transform algorithm developed by Foster (1996). We do not 
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detect any clear signature in the frequency/time spectrum, 
apart from the planet orbital signal. Kepler-10 is intrinsically 
quiet and any stellar activity remains nominal over the course 
of these Kepler observations. The frequency/time spectrum 
does not reveal quarter-dependent fluctuations. 

We then conduct a model-independent Bayesian analysis 
of the entire dataset by employing the Markov Chain Monte 
Carlo (MCMC) implementation presented in Gillon et al. 
(2012). We assume a circular orbit, and set the occultation 
depth, phase-curve amplitude, phase-curve peak offset, pe¬ 
riod, transit duration, time of minimum light and impact pa¬ 
rameter as jump parameters. We assume a simple trapezoid 
function for the occultation, and a Lambertian sphere mod¬ 
ulation for the phase curve. We further assume a quadratic 
law for the limb-darkening (LD) and use ci = 2u\ + U 2 and 
C 2 = u\ - 2 m 2 as jump parameters, where u\ and U 2 are the 
quadratic coefficients. u\ and U 2 are drawn from the theoret¬ 
ical tables of Claret & Bloemen (2011) for the correspond¬ 
ing effective temperature and log g values from Batalha et 
al. (2011). We add the two LD combinations ci and C 2 as 
Gaussian priors in the MCMC fit, based on the theoretical ta¬ 
bles. We run two Markov chains of 10^ steps and assess their 
convergence using the statistical test from Gelman & Rubin 
(1992). We also explore the effects of the data reduction pa¬ 
rameters on deriving the transit parameters. We increase the 
CBV vectors to 8 and reduce and analyze the data separately, 
by pairs of quarters. During all these analyses variations, our 
MCMC fits result in individual transit parameter values within 
1-cr of the final Q1-Q17 values stated above. 

We examine the robustness of the planetary phase-curve 
signal. We stack approximately three years of data, meaning 
that for stellar contamination to happen, stellar activity has to 
be phased exactly on the planetary orbital period of 0.8 days 
(or a multiple). Kepler-10 is an evolved star that is unlikely to 
have vsini consistent with the short orbital period of Kepler- 
10 b. Kepler-10 b is not massive enough to cause any appre¬ 
ciable ellipsoidal or beaming components in the light curve. 
We therefore conclude that the orbital phase-curve is of plan¬ 
etary origin. 

We derive an occultation depth of 7.5 ± 1.5 ppm and a 
phase-curve amplitude of 8.5 ± 1.2 ppm (Figure |5]l. We do 
not detect a phase offset for Kepler-10 b, with the phase off¬ 
set angle constrained to be 9 ± 6 degrees. Our value of the 
occultation depth is 1 - ct consistent with previously reported 
values (5.8±2.5 ppm, Batalha et al. 2011; 9.9± 1.0 ppm, 
Fogtmann-Schulz et al. 2014) 

The occultation depth translates to a brightness temperature 
at the Kepler band of 3220j;[[o K if the planet’s flux is from 
thermal emission, or a geometric albedo of 0.55 ± 0.11 if the 
planet’s flux is from reflection. The phase curve magnitude, if 
solely attributed to reflection, corresponds to an effective geo¬ 
metric albedo of 0.63 ±0.09, or a Bond albedo of 0.94 ±0.13 
for a Lambertian sphere. Given that the zero-albedo dayside- 
average equilibrium temperature of the planet is 2570 K, both 
thermal emission and reflection can have considerable contri¬ 
bution to the planet’s emerging radiation. 

The photometric measurements immediately outside the 
transits contain information about the planet’s nightside. The 
transits occur from -20 to 20 degree (first contact) for this sys¬ 
tem, and we take a 20-degree interval on both sides of the tran¬ 
sit to derive the planet’s nightside emission flux. The result is 
a flux of -0.7 ± 1.2 ppm, which places the 1 - tr upper limit of 
the nightside temperature at 2270 K. Unlike Fogtmann-Schulz 


Total = 

Symmetric Reflection + Asymmetric Reflection +Thermal Emission 



Phase Angle 

Figure 5. Phase-folded light curve of Kepler- 10 b system and examples of 
model fit. Kepler photometry obtained during the quarters 1 to 17 is used to 
derive this light curve, and key steps of data analysis are provided in § 5. 
The light curve is binned to 200 bins for clarity. The magenta lines are 
the modeled phase curves, and the colored lines show contribution of ther¬ 
mal emission (red), symmetric reflection (green), and asymmetric reflection 
(blue). The upper panel shows the simplest model assuming a homogeneous 
surface. The model is dominated by the symmetric reflection component, and 
provides a good fit to the observed phase curve. The lower panel shows the 
best-fit model of the general fit, in which a thermal emission component from 
a hot spot shifted westward has significant contribution. Multiple models can 
provide adequate fit to the phase curve of Kepler- 10 b. 

et al. (2014), our analysis does not yield a definitive nightside 
emission flux. Note that this flux constraint is sensitive to the 
limb-darkening parameterization of the star, as well as the in¬ 
terval size used in the estimate. In any case, the nightside 
brightness temperature must be much lower than the dayside 
brightness temperature. 

5.3. Model-Assisted Analysis 

We then apply our semi-analytical model to analyze the 
phase curve of Kepler-10 b. The eccentricity of the planet’s 
orbit is consistent with zero (Batalha et al. 2011). We per¬ 
form a general fit, in which all parameters are allowed to vary 
in their physically plausible range, and an atmosphere-less fit, 
in which the greenhouse factor is force to / = 1, assuming 
that the planet does not an atmosphere. The atmosphere-less 
fit is motivated by the large bulk density of the planet, and the 
strong irradiation received by the planet (Batalha et al. 2011). 
The posterior distributions resulted from these fits are shown 
in Figure |6] and examples of model fits are shown in Figure 
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Condensation Temperature Cloud Reflection Boosting Factor k 



Heat Redistribution Efficiency e 


General Fit 
Atmosphere-less Fit 


Figure 6. Posterior probability distribution of the parameters for Kepler-10 b from MCMC simulations to fit the phase curve. Black lines show the results of the 
general fit, and the red lines show the results of the atmosphere-less fit in which / = 1 is assumed (i.e., no greenhouse effect). The eclipse depth and the phase 
amplitude derived from the atmosphere-less fit have an upper limit of 9 ppm, and this upper limit corresponds to Ab = 1. This upper limit does not exist when a 
greenhouse effect is allowed in the model. The general fit shows a slight preference to the model has shows a post-eclipse phase shift, driven by the negative heat 
redistribution efficiency. However, this phase offset is not statistically significant. 
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We find that any model with a Bond albedo greater than 0.8 
would provide an adequate fit to the observed phase curve, re¬ 
gardless of whether a greenhouse factor is used, or whether 
asymmetric reflection exists. Figure |5] shows that the best- 
fit model phase curve is dominated by the reflection com¬ 
ponent. The only parameter constrained by the phase curve 
is the Bond albedo, and its value is 0.91;!;o jy, consistent be¬ 
tween the two fits. This value is also consistent with the 
model-independent estimate. This example shows that the 
constraints on the model parameters would be poor for rocky 
planets like Kepler-10 b, because the signal-to-noise ratio of 
the phase amplitude of the planet is only 7-a, much smaller 
than that of Kepler-7 b (30- a). 

Comparing the two fits, the phase curve does not yield ad¬ 
ditional constraints on whether the planet has an atmosphere. 
If Kepler-10 b does not have an atmosphere, its surface must 
be almost fully reflective (i.e., having a reflectivity of almost 
unity) in order to produce the magnitude of the secondary oc- 


cultation and phase curve. Our estimate of the surface albedo 
is consistent with the high end of the range derived by Rouan 
et al. (2011). This is due to the new value of the occultation 
depth reported in this paper, which has a smaller error. The 
linear relationship between the Bond albedo and the occulta¬ 
tion depth in Rouan et al. (2011) is still valid - plugging-in 
our new value of the occultation depth to that equation would 
yield an estimate of the Bond albedo consistent with ours. 
Such a high reflectivity makes a solid surface possible. The 
day side average equilibrium temperature would be 1720 K for 
a Bond albedo of 0.8, and only 1450 K for a Bond albedo of 
0.9, and these temperatures are lower than the melting temper¬ 
ature of some materials such as aluminum oxide. Therefore 
it is unclear whether Kepler-10 b has a molten lava surface as 
suggested by Rouan et al. (2011); but if it does, the molten 
lava has to be highly reflective. 

If Kepler-10 b has a thin atmosphere that provides moder¬ 
ate greenhouse effect (/ > 1), the thermal emission compo¬ 
nent can contribute to the planetary flux significantly. We And 
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that the fit to the observed phase curve can be slightly im¬ 
proved, if a hot spot shifted to the west side of the substellar 
point is included (Figure |5]l. In the best-fit model, the shifted 
hot spot creates a post-eclipse phase offset that appears to be 
more consistent with the observations. The posterior distribu¬ 
tions also show a preference to this class of solutions (Figure 
However, due to the low signal-to-noise ratio of the phase 
offset, this improvement is marginal, and the constraint on the 
greenhouse factor is poor. 

6. DISCUSSION 

6.1. Two Solutions for a Phase Ojfset 

When applying our semi-analytical model to study the 
phase curve of Kepler-7 b, we find two classes of solution 
can fit the phase curve equally well. The two classes of solu¬ 
tion cannot be distinguish via usual statistical diagnostics like 
or the Bayesian Information Criterion. Here we generalize 
this dual solutions to any visible-wavelength phase curve that 
shows a significant phase offset for short-period exoplanets, 
and discuss how to break the degeneracy. 

Take the post-eclipse phase offset for an example. Such an 
offset, qualitatively, means that the planet’s brightness (from 
thermal emission, reflection, or both) of the west side of the 
substellar point is greater than that of the east side. If the 
phase curve is principally produced by reflection, the west 
side of the substellar point on the planet would have a greater 
albedo than the east side. If the phase curve is principally 
produced by thermal emission, the west side of the substellar 
point would be hotter than the east side. Either due to reflec¬ 
tion or thermal emission, the phase shift indicates east-west 
transport on the planet. Here we see that, although the phase 
shift implies transport, there can be ambiguity in interpreting 
the direction of the phase shift, mainly due to entanglement 
of thermal emission and reflection. 

Atmospheric circulation models may provide further con¬ 
straints to break the degeneracy. For tidally-locked irradiated 
planets having a wide range of temperatures, atmospheric cir¬ 
culation models uniformly predict super-rotating equatorial 
winds and eastward shifted hot spot (e.g.. Showman & Guil- 
lot 2002; Showman et al. 2011). For these planets, a pre¬ 
eclipse phase offset would be driven by thermal emission, and 
the atmosphere can be homogeneously reflective; and a post¬ 
eclipse phase offset would be driven by asymmetric reflection 
due to formation of patchy clouds. 

Another way to break the degeneracy is to measure the 
planet’s flux at a different wavelength, using different wave¬ 
length dependencies between the reflection and the thermal 
emission component. In the example of Kepler-10 b, a sin¬ 
gle measurement of the secondary occultation depth at mid- 
infrared wavelengths can rule out the scenarios in which ther¬ 
mal emission drives the phase offset. One caveat is that the 
opacities of the atmosphere can be very different between 
the visible and longer wavelengths, and the thermal emission 
at these wavelengths may probe the temperature at different 
pressure levels. Using planetary atmosphere chemistry and 
radiative transfer models can potentially mitigate this caveat 
(e.g., Hu & Seager 2014). The joint observations at infrared 
wavelengths can only be done for few Kepler gaseous plan¬ 
ets, due to the faintness of the sources. Looking ahead, TESS 
and PLATO will measure visible-wavelength phase curves of 
gaseous and rocky planets around nearby stars. Eor those sys¬ 
tems, the James Webb Space Telescope can obtain their sec¬ 
ondary occultation depths and provide great constraints to dis¬ 


tinguish the two classes of models, and pinpoint the albedo 
and the thermal energy transfer of the planets. 

Once this degeneracy can be broken with additional mod¬ 
eling or observational efforts, a lot more can be known about 
the planets. Eor example, if the phase offset is principally pro¬ 
duced by thermal emission, the greenhouse factor / could be 
tightly constrained by the phase curve. Eor a terrestrial exo¬ 
planet that may or may not have an atmosphere, a significant 
constraint of / > 1 would indicate the planet must have an 
atmosphere. Eor another example, if the phase offset is prin¬ 
cipally produced by reflection, the cloud condensation tem¬ 
perature or the surface freezing temperature could be derived, 
which would indicate cloud or surface materials. 

In this work we have assumed that all surface elements are 
Lambertian surfaces. For homogeneous scenarios, this as¬ 
sumption leads to Ag = 2Ag/3. The real planetary atmospheres 
may be significantly backscattering. For example, if a planet 
has a semi-infinite atmosphere that is purely Rayleigh scatter¬ 
ing, Ag = 3Ab/ 4 (Dlugach & Yanovitskij 1974). By neglect¬ 
ing backscattering we might overestimate the Bond albedo, 
by 10- 15%. Correcting this bias leads to 1% decrease in 
the estimate of the greenhouse factor / for the homogeneous 
atmosphere solutions, or 1% increase in the estimate of the 
cloud condensation temperature Tc for the patchy cloud solu¬ 
tions. 

In this paper we focus on planets that have circular orbits, 
which would include most of the planets whose phase curves 
can be measured by Kepler (Esteves et al. 2014). Eor com¬ 
pleteness, we note that planets that have eccentric orbits may 
have phase offsets without heat redistribution (Cowan & Agol 
2011; Kataria et al. 2013). With heat redistribution, the phase 
curves of eccentric planets may show more complex patterns 
than the planets in circular orbits. 

6.2. Formation of Patchy Clouds 

The patchy cloud scenario is the leading explanation for the 
post-eclipse phase offset of Kepler-7 b. The inferred pattern of 
clouds follows a longitudinal distribution controlled by the lo¬ 
cal temperature of the atmosphere. AnnendixlBlshows that the 
condensation and sublimation are the dominant microphysical 
processes for 1-10 micron-sized particles at a temperature of 
1000-2000 K and a pressure of 0.1-1 bars. Under these condi¬ 
tions, advection cannot transport cloud particles fast enough 
- compared with evaporation - to All the “hole” in the cloud. 
Therefore, the patchy clouds may be condensational clouds 
that form instantaneously when the temperature falls below 
the condensation temperature. 

In this picture, the thermal hot spot at 0.1-1 bars is the con¬ 
trolling dynamical pattern that affects the cloud appearance 
and the phase curve. However, at lower pressures, e.g. 1 
mbar, the circulation would be dominated by symmetric flows 
from the substellar point to the anti-stellar point (e.g.. Show¬ 
man et al. 2008, 2009). As a result, the 0.1-1 bar thermal hot 
spot may correspond to a low temperature at 1 mbar, which 
could lead to another cold trap. The question is, can this low- 
pressure cold trap lead to cloud formation? We may consider 
this question by comparing the two timescales: the advective 
timescale of 10"^ s (Appendix iBl). and the vertical diffusive 
timescale. The latter can be estimated as T/^/Kzz, where H 
is the pressure scale height and Kzz is the vertical diffusiv- 
ity. Using H = 50 km and Kzz = 5 x 10* cm^ s“' (Parmentier 
et al. 2013), we estimate the vertical diffusive timescale to 
be 5 X (O'* s. Therefore, on these planets, the condensable 
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species can only be transported upward by a fraction of the 
pressure scale height for a typical advection timescale. The 
picture of the patchy clouds can then be, that on the nightside 
and the west side of the substellar point, the planet has a cold 
trap for the condensable species at the pressure of 0.1-1 bar, 
and due to the cold trap, the species is depleted at lower pres¬ 
sures. On the east side of the substellar point, the species can 
be transported upward, but this vertical transport would be ef¬ 
fectively quenched by fast zonal winds and then clouds would 
not form at the low-pressure cold trap. Such a scenario can be 
studied by a model coupling atmospheric dynamics and cloud 
microphysics that includes condensation and sublimation. 

The conditions to form the patchy clouds would thus be the 
availability of condensable material, a suitable atmospheric 
temperature, and a fast zonal wind. These conditions are 
met on gaseous planets having an equilibrium temperature 
of 1000-2000 K. As a result, patchy clouds featuring strong 
longitudinal variation of cloud coverage may be a common 
phenomena in the atmosphere of irradiated gaseous exoplan¬ 
ets. Although the cloud condensation temperature is fixed 
for a certain material, a wide range of stellar irradiation can 
still lead to patchy clouds, because the clouds can form at 
appropriate pressures corresponding to the condensation tem¬ 
perature. Our semi-analytical model has also been applied 
to two other Kepler planetary systems that show phase varia¬ 
tion which cannot be due to orbital effects (i.e., doppler beam¬ 
ing and ellipsoidal effects). Similar post-eclipse phase offsets 
due to displaced patchy clouds are found, but diversities in 
the atmosphere and cloud reflectivities emerge (Shporer & Hu 
2015). 

7. CONCLUSION 

To analyze the visible-wavelength phase curves of exo¬ 
planets, we develop a semi-analytic model to calculate the 
reflection and thermal emission components of exoplanets 
as a function of orbital phase. The phase curve model de¬ 
pends on five parameters that characterize key energetic pro¬ 
cesses of the planets: the Bond albedo, the heat redistri¬ 
bution efficiency, the greenhouse effect factor (set to unity 
for atmosphere-less models), the condensation temperature 
(not used for symmetric reflection models), and the reflec¬ 
tivity boosting factor (not used for symmetric reflection mod¬ 
els). The model tremendously simplifies the otherwise time- 
consuming computation for atmospheric circulation, cloud 
formation, and radiative transfer, which enables exploration 
of the full parameter space by the MCMC method. 

Our model framework applied to Kepler-7 b reveals a gen¬ 
eral degeneracy in interpreting the single-band phase curves 


of a hot exoplanet. If both reflection and thermal emission 
contributes to the planetary flux in the wavelength of interest, 
a potential phase shift can be attributed to a hot spot shifted 
to one direction in a clear atmosphere, or a hot spot shifted to 
the opposite direction in an atmosphere that has patchy clouds 
or on a surface that has inhomogeneous reflectivities. For 
Kepler-7 b we find both scenarios fit the visible-wavelength 
phase curve fairly well, despite vastly different atmospheric 
circulation regimes they imply. For Kepler-7 b, atmospheric 
circulation models and previous Spitzer observations rule out 
the model attributing the phase offset to the thermal emission 
component, and make the patchy cloud scenario the only plau¬ 
sible explanation. Our model can therefore derive strong con¬ 
straints on the reflectivities of the clear part and the cloudy 
part of the atmosphere, and the location of the cloud bound¬ 
ary. We also present a full phase curve of Kepler- 10 b, and im¬ 
proved measurements of the occultation depth and the phase 
magnitude. We do not detect a significant phase offset in the 
phase curve of Kepler-10 b. Due to the low signal-to-noise 
ratio of its phase amplitude, our model can only constrain the 
Bond albedo of the planet. In all, the phase curve analysis is 
a powerful tool to investigate the exoplanets discovered by 
Kepler and those to be discovered by future transit survey 
missions. This “cost-free” information from current planet 
searching strategy will yield highly valuable constraints on 
the planets’ characteristics. 
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APPENDIX 

PHASE FUNCTION ASSOCIATED WITH ASYMMETRIC DISTRIBUTION OF CONDENSED PARTICLES 

The local longitude range on the planet illuminated by the star is [“ f, ^] ■ For the scenarios with asymmetric reflection, the 
local longitudes ranging in [^ 1 ,^ 2 ] have a lower reflectivity than other longitudes on the planet. and ^2 are set by comparison 
between the local temperature and the freezing temperature of the surface or the condensation temperature of the condensable 
species in the atmosphere, and they must also obey - f <6<6<f- 
The phase function of the reflected light associated with the additional reflection, which is in addition to the reflected light by 
a homogeneous surface or atmosphere, is 



TT 


cos(a-(p) COS (j)d(j), 


/ solid 


(Al) 


in which the integration covers the longitudes that have the high reflectivity and that can be seen along the line of sight. Given 
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the relationship between the local longitude and the observer longitude, ^ = (j) — a, the phase function can also be integrated as 


4)= — / cos^cos(^ + a)d^. 


(A2) 


/ solid 


Define a' = -a, and Equation (IA2t with respect to a' takes the same form as Equation (lAll) with respect to a. It is convenient to 
evaluate the phase function in a' and the exact formulae for this integral for the phase function are given as follows. 

Eor -TT < a' < 0, 

If-f < a' + f <6, 


$ = — [cos a'in + a') + sin(7r + a')] ; 


(A3) 


If6 <a' + f <6, 


If6<a' + f <f 


Eor TT > a' > 0, 

If-f < a'-f <6, 


1 TT 1 

$ = — {cosa'(—+ ^i) + -[sin(7r+a') + sin(2^i - a')]}; 
TT 2 2 


4) = — [cos a'(7r + a' + - ^ 2 ) + sin(7r + a') + cos(^ i+^i-cx') sin(^ 1 - ■C 2 )] i 

TT 


$ = — [cosa'(7r-Q;' + ^i -^2) + sina^ + cos(^i + ^2-Q!0sin(^i -^ 2 )]; 


(A4) 


(A5) 


(A6) 


if6 <a'-f <a 


If6<a'-f <f, 


1 TT 1 

= - {cos a'( - - 6) + T [sin a' - sin( 2^2 - a ')]}; 

TT A A 


$ = —[cosQ!'(7r-Q!0 + sina']. 

TT 


(A7) 


(A8) 


We also need to calculate the contribution of the additional asymmetric reflection to the planet’s Bond albedo. This is done by 
further integrating the phase function with respect to a. The phase integral is 


?'(6,6)= [ $sin|Q!|^fa. 

J-TT 


(A9) 


The value of q' as a function of and ^2 is evaluated numerically and shown in Eigure0 Eor example, if cloud particles cover 
the whole substellar hemisphere (^2 = Ci)’ 4 - 3/2, which gives the classical result of a Lambertian sphere. The Bond albedo of 
the planet is 

Ab = rQ+-q (AlO) 



Figure 7. The phase integral associated with the frozen surface or the condensed particles {q') as a function of and ^ 2 - 
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MICROPHYSICAL CONSIDERATION OF THE PATCHY CLOUDS 

To simplify the lava lake scenario and the patchy cloud scenario, we assume that the pattern for lava or clouds follows a longi¬ 
tudinal distribution controlled by the local temperature of the surface or the atmosphere. For the patchy clouds, this is equivalent 
to saying that condensation and sublimation are the dominant microphysical processes, while precipitation and advection are 
relatively unimportant. Here we show that this assumption is likely to be valid for 1-10 micron-sized particles at a pressure level 
of 0.1 bars with a temperature of ^ 1000 K. With this assumption, a non-uniform longitudinal distribution of temperature, as a 
result of energy redistribution, can lead to a non-uniform longitudinal distribution of cloud particles. 

For the patchy cloud scenario we have assumed that condensation and sublimation is controlled solely by the local temperature. 
To investigate when this “local cloud” assumption is valid, we compare the following timescales: 

• Advective timescale Tadv ~ Rp/u, where Rp is the radius of the planet and u is the equatorial wind speed, u would be on 
the order of 1 km s“' for irradiated exoplanet atmospheres (e.g.. Showman & Guillot 2002). For Kepler-10 b Tadv ~ 10"^ s, 
and for Kepler-7 b it would be longer. 

• Condensation/sublimation timescale Tcond i 2 gnf(p-p.i) ’ where Dp is the particle diameter, pp is the density of the condensed 
material, kb is the Boltzmann constant, T is the dayside temperature, K is the diffusivity of the atmosphere, m is the 
molecular mass of the condensable species, and p - psat is the difference between the partial pressure and the saturation 
vapor pressure of the condensable material. The estimation formula is derived from the condensational growth equation 
of atmospheric aerosols for perfect accommodation and the continuum regim^^ (Seinfeld & Pandis 2006). Note that the 
condensation timescale is a strong function of the particle size and the partial pressure of the condensable species. Estimates 
of the condensation timescale using typical values are shown in Table 2. 


Table 2 

Condensation timescale of a potentially condensable species in the atmosphere. The typical values are taken as pp = 3 X 10^ kg m“^ (enstatite), T = 2000 K, 

D = 6x 10“^ m- s“*, m = L7 X 10“^^ kg, and atmosphere pressure P = 10^ Pa. 



Dp = m 

Dp = 10^5 m 

p-ftat = 10-^ P 

70s 

7 X 10^ s 

P-Psat = P 

7 X 10^ s 

7 X 10^ s 


• Gravitationally settling timescale Tjettie 7 ^^ where H is the atmospheric scale height and Vsettie is the settling velocity. 
At 10"^ Pa, the settling velocity for a 1 pm sized particle is in the order of 0.01 m s“', and the atmospheric scale height on 
Kepler-10 b is about 30-80 km. We estimate Tsettie ~ 3-8 x 10^ s. This value is greater by a factor of a few for Kepler-7 b. 
For 10 pm sized particle Tsettie would be smaller by 2 orders of magnitude. 

Based on these timescale estimates, we find it possible to have Tcond < Tadv < Tsettie for particles of 1-10 micron radius at 
the pressure level of 0.1 bars in the atmospheres of Kepler-10 b and Kepler-7 b. That is to say, formation and sublimation 
of a potential condensable species is in local equilibrium, i.e., the amount of condensate particles only depends on the local 
temperature. In comparison, horizontal and vertical transport can be much less efficient. These estimates also suggest that the 
particle size could be greatly limited by advection in the atmospheres having strong zonal winds, instead of limited by vertical 
transport as for some brown dwarfs (e.g.. Cooper et al. 2003). 
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